En primer lugar se cargan los paquetes necesarios para el análisis, y se carga la tabla de datos original “Datos de uso Bicimad Agosto 2018”.
# Load de paquetes para analisis
library(grDevices)
library(geosphere)
library(tidyverse)
library(data.table)
library(plotly)
library(factoextra)
library(ggpubr)
library(Hmisc)
library(corrplot)
library(cluster)
library(NbClust)
library(reprex)
library(magrittr)
library(readr)
library(Rtsne)
# Load dataset target del estudio
# Datos originales de Agosto 2018. 274 mil registros
data_orig <- fread('Bicimad_201808_3.csv', na.strings = '\\N')Utilizamos datos externos para enriquezer el dataset. Se han utilizado las siguientes fuentes:
Tb. Estaciones - Dataset que contiene el id de estación, dirección, código postal, distrito, coordenadas geográficas… Estos datos se han enriquecido para las estaciones tanto de origen como destino.
Tb. Meteorológica - Datos metereológicos de Madrid por fecha con datos de temperatura y rachas de viento.
Tb. Edad - Tabla que relaciona el código de edad con la descripción literal.
Tb. Tipo de Usuario - Tabla que relaciona el código de tipo de ususario con su descripción literal.
##### LECTURA DE TABLAS #####
# Tb. estaciones -- inicio. Se carga el maestro de estaciones.
master_estaciones_ini <- fread('tb_stations.csv', colClasses = 'character')
master_estaciones_ini <- master_estaciones_ini %>%
rename(name_station_ini = name,
longitude_station_ini = longitude,
latitude_station_ini = latitude,
address_station_ini = address,
idunplug_station_ini = id,
distrito_station_ini = distrito,
cp_ini = cod_postal) %>%
dplyr::select(-number)
# Tb. estaciones -- fin.
master_estaciones_fin <- fread('tb_stations.csv', colClasses = 'character')
master_estaciones_fin <- master_estaciones_fin %>%
rename(name_station_fin = name,
longitude_station_fin = longitude,
latitude_station_fin = latitude,
address_station_fin = address,
idplug_station_fin = id,
distrito_station_fin = distrito,
cp_fin = cod_postal) %>%
dplyr::select(-c(number))
# Tb. Meteo. Se carga y mapea el dataset con los datos climatológicos de agosto.
master_meteo <- fread('tb_meteo.csv')
master_meteo <- master_meteo %>%
rename(unplug_date = fecha) %>%
mutate(temp_max = as.numeric(sub(",", ".", temp_max, fixed = TRUE)),
temp_max_hora = as.numeric(sub(",", ".", temp_max_hora, fixed = TRUE)),
temp_min = as.numeric(sub(",", ".", temp_min, fixed = TRUE)),
temp_min_hora = as.numeric(sub(",", ".", temp_min_hora, fixed = TRUE)),
temp_media = as.numeric(sub(",", ".", temp_media, fixed = TRUE)))
# Tb. Edad. Maestro de edad
master_edad <- fread('tb_ageRange.csv', colClasses = 'character')
# Tb. userType. Maestro de tipo de usuarios
master_user <- fread('tb_userType.csv', colClasses = 'character')En el dataset enriquecido se ha filtrado por código de estación para trabajar únicamente con las estaciones en el distrito “Distrito Retiro”, añadiendo estaciones que por proximidad pueden resultar interesante spara el estudio.
El dataset enriquezido, de ahora en adelante referido como dataset target (var(data)), se ha cruzado por campos comunes y se han normalizados los data types. Además, se ha calculado la distancia mínima entre estaciones computando la Haversine distance entre coordenadas.
# DATA WRANGLING
# Filtrar por estaciones que nos interesan
estaciones_retiro <- c(64,65,66,69,72,73,74,75,78,79,80,84,85,86,88,90,91,99,100,101,102,107)
# Se unen dataset con los maestros de estaciones, datos meteorológicos, edad y tipo de usuario
data <- data_orig %>%
rename(movimiento_id = "id_mov") %>%
filter(idunplug_station_ini %in% estaciones_retiro | idplug_station_fin %in% estaciones_retiro) %>%
mutate(idunplug_station_ini = as.character(idunplug_station_ini),
cp_ini = as.character(cp_ini),
idunplug_base = as.character(idunplug_base),
idplug_station_fin = as.character(idplug_station_fin),
cp_fin = as.character(cp_fin),
idplug_base = as.character(idplug_base),
user_type = as.character(user_type),
user_age_range = as.character(user_age_range),
user_zip_code = as.character(user_zip_code),
bank_holiday_flg = as.character(bank_holiday_flg),
unplug_hour = as.character(unplug_hour),
unplug_date=as.character(unplug_date)
) %>%
dplyr::select(-c(unplug_hourtime)) %>%
left_join(master_edad, by='user_age_range') %>%
left_join(master_user, by='user_type') %>%
left_join(master_estaciones_ini, by=c('idunplug_station_ini', 'cp_ini')) %>%
left_join(master_estaciones_fin, by=c('idplug_station_fin', 'cp_fin')) %>%
mutate(longitude_station_ini = as.numeric(longitude_station_ini),
latitude_station_ini = as.numeric(latitude_station_ini),
longitude_station_fin = as.numeric(longitude_station_fin),
latitude_station_fin = as.numeric(latitude_station_fin)) %>%
mutate(distance = distHaversine(cbind(longitude_station_ini, latitude_station_ini), cbind(longitude_station_fin, latitude_station_fin))) %>%
left_join(master_meteo, by='unplug_date')
rm(data_orig, master_edad, master_estaciones_fin, master_estaciones_ini, master_meteo, master_user, estaciones_retiro)
# Ordenar variable 'day'
data$day <- factor(data$day, levels=c('L','M','X','J','V','S','D'))
# save(data, file='data.RData')El posterior análisis exploratorio describe el dataset target por variables individuales y cruzes entre variables mediante análisis gráficos y sus estadíticos principales.
Más allá de que el distrito 02 - retiro sea el mayoritario (pues estamos influyendo en ello), se observa que los distritos más recurrentes como inicio de viaje son los colindantes al propio retiro, aunque destaca que distrito se situe por encima de Salamanca, cuando es este último el distrito más cercano al propio parque, por lo que ya denota una tendencia natural de los usuarios.
x <- data %>% group_by(distrito_station_ini) %>% summarise(n = length(user_id))
rm(x)
x <- data %>%
group_by(distrito_station_ini) %>%
summarise(n = length(distrito_station_ini)) %>%
ungroup() %>%
mutate(porc = n/sum(n)) %>%
filter(porc*100 >= 0.1)
z <- ggplot(x, aes(x=distrito_station_ini, y=porc*100))+theme_classic()+
geom_bar(stat = "identity",fill="lightblue")+ labs(title="Distrito de partida", x="", y="%")+theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplotly(z)Los códigos postales de las estaciones de inicio directamente vinculados con la situación de las estaciones siguen la misma distribución.
x <- data %>%
group_by(cp_ini) %>%
summarise(n = length(cp_ini)) %>%
ungroup() %>%
mutate(porc = n/sum(n))
# p <- ggplot(x, aes(x=cp_ini, y=porc*100)) +
# geom_bar(stat = "identity") +
# theme(axis.text.x = element_text(angle = 45, hjust = 1))
p <- ggplot(x, aes(x=cp_ini, y=porc*100))+theme_classic()+
geom_bar(stat = "identity",fill="lightblue")+ labs(title="Código postal de inicio", x="", y="%")+ theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplotly(p)# 1.2 ID UNPLUG_STATION INI
x <- data %>% group_by(idunplug_station_ini) %>% summarise(n = length(user_id))
x <- data %>%
group_by(idunplug_station_ini) %>%
summarise(n = length(idunplug_station_ini)) %>%
ungroup() %>%
mutate(porc = n/sum(n)) %>%
filter(porc*100 >= 0.1)
p <- ggplot(x, aes(x=idunplug_station_ini, y=porc*100))+theme_classic()+
geom_bar(stat = "identity",fill="lightblue")+ labs(title="Estación de partida", x="", y="%") + theme(axis.text.x = element_text(angle = 90))
ggplotly(p)Para los distritos en los que finalizan los trayectos, vuelve a aparecer una tendencia en dirección a la zona centro superando a Salamanca que es la distrito más cercano al parque del retiro.
x <- data %>% group_by(distrito_station_fin) %>% summarise(n = length(user_id))
rm(x)
x <- data %>%
group_by(distrito_station_fin) %>%
summarise(n = length(distrito_station_fin)) %>%
ungroup() %>%
mutate(porc = n/sum(n)) %>%
filter(porc*100 >= 0.1)
z <- ggplot(x, aes(x=distrito_station_fin, y=porc*100))+theme_classic()+
geom_bar(stat = "identity",fill="lightblue")+ labs(title="Distrito de llegada", x="", y="%")+ theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplotly(z)Distribución calcada con las estaciones de inicio, pues estamos eligiendo aquellos viajes que comienzan o terminan en el retiro.
x <- data %>%
group_by(cp_fin) %>%
summarise(n = length(cp_fin)) %>%
ungroup() %>%
mutate(porc = n/sum(n))
# p <- ggplot(x, aes(x=cp_fin, y=porc*100)) +
# geom_bar(stat = "identity") +
# theme(axis.text.x = element_text(angle = 45, hjust = 1))
p <- ggplot(x, aes(x=cp_fin, y=porc*100))+theme_classic()+
geom_bar(stat = "identity",fill="lightblue")+ labs(title="Codigo postal de fin", x="", y="%")
ggplotly(p)Conclusiones similares a las estaciones de inicio, donde al estar filtrando por los viajes que se inician o terminan en el retiro estamos forzando que estas estaciones sean más frecuentes.
x <- data %>%
group_by(idplug_station_fin) %>%
summarise(n = length(idplug_station_fin)) %>%
ungroup() %>%
mutate(porc = n/sum(n)) %>%
filter(porc*100 >= 0.1)
# p <- ggplot(x, aes(x=idplug_station_fin, y=porc*100)) +
# geom_bar(stat = "identity") +
# theme(axis.text.x = element_text(angle = 45, hjust = 1))
p <- ggplot(x, aes(x=idplug_station_fin, y=porc*100))+theme_classic()+
geom_bar(stat = "identity",fill="lightblue")+ labs(title="Estación de llegada", x="", y="%")
ggplotly(p)Los usuarios conabono anual son muy claramente mayoritarios, estando presentes en el 86% de los viajes que se inician o terminan en el retiro, con menos de un 10% los trabajadores de Bicimad y con el 4% los usuarios con abono ocasional.
x <- data %>%
group_by(user_type_lit) %>%
summarise(n = length(user_type_lit)) %>%
ungroup() %>%
mutate(porc = n/sum(n))
# p <- ggplot(x, aes(x=user_type, y=porc*100)) +
# geom_bar(stat = "identity") +
# theme(axis.text.x = element_text(angle = 45, hjust = 1))
p <- ggplot(x, aes(x=user_type_lit, y=porc*100))+theme_classic()+
geom_bar(stat = "identity",fill="lightblue")+ labs(title="Tipo de usuario", x="", y="%")
ggplotly(p)Anque la mayoría de usuarios no declaran la edad, de aquellos que la declaran el rango. Si excluimos a estos y recalculamos los procentajes estaríamos diciendo que la mitad de los viajes con origen o destino en el retiro son por usuarios de entre 27 y 40 años, donde si unimo los usuarios de 27 a 40, tendríamos un peso del 86% de usuarios adultos. Con esto coge peso la teoría de que es un servicio usado por gente adulta, pudiendo ser usada para desplazamientos en dirección al trabajo, por ejemplo.
x <- data %>%
group_by(user_age_range_lit) %>%
summarise(n = length(user_age_range_lit)) %>%
ungroup() %>%
mutate(porc = n/sum(n))
p <- ggplot(x, aes(x=user_age_range_lit, y=porc*100)) +
theme_classic() +
geom_bar(stat = "identity",fill="lightblue") +
labs(title="Rango de edad de usuario", x="Grupo", y="%")+
theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplotly(p)Aunque gran parte de los usuarios no declaran su codigo portal, en el histograma destaca la perdida de peso del código postal 28014, que si estaba muy presente en las estaciones de inicio y fin. Esto podría ser acausa de que esta localización es más turistica que residencial.
x <- filter(data,!is.na(user_zip_code))
x <- x %>%
group_by(user_zip_code) %>%
summarise(n = length(user_zip_code)) %>%
ungroup() %>%
mutate(porc = n/sum(n))%>%
filter(porc*100 >= 0.1)
p <- ggplot(x, aes(x=user_zip_code, y=porc*100)) +
theme_classic() +
geom_bar(stat = "identity",fill="lightblue") +
labs(title="Codigo postal de usuario", x="", y="%")+
theme(axis.text.x = element_text(angle = 90, hjust = 1))
ggplotly(p)En base al histograma, el servicio de bicimad para los viajes relativos a retiro, parece ser usado como media de transporte rutinaria, más que a modo ocasional, con una tendencia muy marcada durante la semana, que cae en el fin de semana. Es posible que si cogiesemos un dataset de otro mes que no sea de un periodo vacacional, estas diferencias sean más marcadas todabía.
Se observa que el 77% de los trayectos ocurren en dias laborables, y el restate en festivos(incluyendo el dia 15 de Agosto).
x <- data %>%
group_by(bank_holiday_flg) %>%
summarise(n = length(bank_holiday_flg)) %>%
ungroup() %>%
mutate(porc = n/sum(n))
p <- ggplot(x, aes(x=bank_holiday_flg, y=porc*100)) +
theme_classic() +
geom_bar(stat = "identity",fill="lightblue") +
labs(title=" Laborables vs Festivos", x="", y="%")
ggplotly(p)Se notan horas pico coincidiendo con el inicio (8) y fin (19) de la que podría ser una jornada laboral, acentuando además las 15h entendiendolo como la hora de comer. En un análisis más profundo, filtrando por bank_holiday_flg = 0, es decir, dias laborables, esta diferencia se ve acentuada.
x <- data %>%
group_by(unplug_hour,bank_holiday_flg) %>%
summarise(n = length(unplug_hour)) %>%
ungroup() %>%
mutate(porc = n/sum(n))
x$unplug_hour=as.numeric(x$unplug_hour)
p <- ggplot(x, aes(x=unplug_hour, y=porc*100)) +
theme_classic() +
geom_bar(stat = "identity",fill="lightblue") +
labs(title="Hora de inicio del trayecto", x="", y="%") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggplotly(p)Transformando la variale numérica de travel_time a minutos para hacerla más legible, observamos que aunque existen outliers, con viajes de 8.000 min o de 0.17 min, la media de los viajes está entorno a 11 minutos, situándose el 3er cuartil por debajo de los 20 minutos.
Descripcion de estadisticos:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.017 7.817 11.400 20.827 18.050 8030.233
Desviacion estandar:
## [1] 87.87702
Varianza:
## [1] 7722.371
Observamos que la distancia media recorrida está entorno a los 2KM. La mayoría de los viajes están contenidos en la franja de los 1-3km.
Descripcion de estadisticos:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 1235 1823 1889 2397 5944
Desviacion tipica:
## [1] 956.4534
Varianza:
## [1] 914803.1
Altas temperaturas sin bajas de los 30 grados, sin haber grandes dispersiones. De haber sido otros meses para analizar, habría cogido valor, el introducir datos referentes a las precipitaciones
Descripcion de estadisticos:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 30.40 32.70 33.80 34.43 36.70 40.00
Desviacion tipica:
## [1] 2.508407
Varianza:
## [1] 6.292106
De igual manera las temperaturas minimas con poca dispersión e igualmente altas.
Descripcion de estadisticos:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 16.80 19.60 21.00 21.36 23.00 25.90
Desviacion tipica:
## [1] 2.341467
Varianza:
## [1] 5.48247
Descripcion de estadisticos:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 11.00 13.00 15.00 16.26 19.00 25.00
Desviacion tipica:
## [1] 3.597168
Varianza:
## [1] 12.93962
En primera instancia, cuando tratamos de analizar la variable travel time, los outliers no nos dejan ver realmente como se comportan los usuarios, en el momento que los dejamos de lado para pintar el boxplot, observamos lo que nos desprendió el análisis estadistico, donde el servicio era usado para viajes relativamente cortos en su mayoría.
# 1.C.1 TRAVEL_TIME
# Anadir 'travel_time' en minutos
data <- data %>%
mutate(travel_time_mins = travel_time/60)
# Identificacion de outliers
outliers_travel_time_mins <- boxplot.stats(data$travel_time_mins)$out
p <- ggplot(data, aes(y=travel_time_mins)) + geom_boxplot() + labs(title="Travel_time - outliers")
ggplotly(p)# Quitando outliers de 'data'
data_no_outliers <- data[-which(data$travel_time_mins %in% outliers_travel_time_mins),]
# save(data_no_outliers, file='data_no_outliers.RData')
# Boxplot - 'data_no_outliers'
p <- ggplot(data_no_outliers, aes(y=travel_time_mins)) + geom_boxplot() + labs(title="Travel_time - no outliers")
ggplotly(p)Los outliers no dificultan tanto el análisis de la variable en este caso, aun así en el momento que los excluimos observamos que en terminos generales los viajes rondan los 1 a 3 KM.
# 1.C.2 DISTANCE
# Identificacion de outliers
outliers_distance <- boxplot.stats(data$distance)$out
p <- ggplot(data, aes(y=distance)) + geom_boxplot() + labs(title="Distance - outliers")
ggplotly(p)Para conocer la relación entre dos variables categóricas utilizaremos:
Tablas de contingencia para describir esta relación a nivel muestral. Test de chi-cuadrado para hacer inferencia sobre la relación entre ambas variables.
#0 - Unknown
#1 - 0 y 16 años
#2 - 17 y 18 años
#3 - 19 y 26 años
#4 - 27 y 40 años
#5 - 41 y 65 años
#6 - > 66 años
frec <- round(prop.table(table(data$user_type, data$user_age_range))*100, 2)
frec##
## 0 1 2 3 4 5 6
## 1 34.83 0.20 0.71 5.56 27.53 17.37 0.50
## 2 3.83 0.00 0.00 0.02 0.01 0.01 0.00
## 3 3.49 0.74 0.00 0.28 2.34 2.59 0.00
##
## Pearson's Chi-squared test
##
## data: table(data$user_type, data$user_age_range)
## X-squared = 6806.9, df = 12, p-value < 2.2e-16
#Visualizamos la relación en un histograma
#heatmap(table(data$user_type, data$user_age_range), Colv = NA, Rowv = NA, xlab="user_age_range", ylab="user_type")
barplot(frec, beside=FALSE,legend.text= rownames(frec[1]), main='Distribución % Tipo usuario por edad', xlab='Rangos Edad', ylab='% Distribución tipo usuario')#Generamos la tabla de variables eliminando como distrito inicial el Retiro
table_distrit_age <- data %>%
dplyr::select(c(distrito_station_ini,user_age_range)) %>%
filter(grepl("03",distrito_station_ini)==FALSE) %>%
mutate(distrito_station_ini = as.factor(distrito_station_ini),
user_age_range = as.factor(user_age_range))
#Realizamos la matriz de contingencia de las frecuencias relativas
round(prop.table(table(table_distrit_age))*100, 2)## user_age_range
## distrito_station_ini 0 1 2 3 4 5 6
## 01 CENTRO 17.04 0.21 0.25 2.17 11.90 6.68 0.10
## 02 ARGANZUELA 7.56 0.25 0.13 1.06 4.93 3.28 0.04
## 04 SALAMANCA 13.01 0.26 0.19 1.64 9.57 6.05 0.21
## 05 CHAMARTÍN 1.01 0.04 0.05 0.24 0.88 0.69 0.00
## 06 TETUÁN 0.52 0.02 0.00 0.06 0.39 0.24 0.00
## 07 CHAMBERÍ 3.00 0.04 0.04 0.47 2.11 1.59 0.03
## 09 MONCLOA 0.86 0.00 0.02 0.08 0.65 0.43 0.01
#Incluimos simulación Chi-cuadrado
chisq.test(table_distrit_age$distrito_station_ini,table_distrit_age$user_age_range)##
## Pearson's Chi-squared test
##
## data: table_distrit_age$distrito_station_ini and table_distrit_age$user_age_range
## X-squared = 215.74, df = 36, p-value < 2.2e-16
#Obtenemos un Warning al tener valores inferiores a 5 y realizamos también el Test de fisher
fisher.test(table_distrit_age$distrito_station_ini,table_distrit_age$user_age_range,simulate.p.value=TRUE)##
## Fisher's Exact Test for Count Data with simulated p-value (based
## on 2000 replicates)
##
## data: table_distrit_age$distrito_station_ini and table_distrit_age$user_age_range
## p-value = 0.0004998
## alternative hypothesis: two.sided
#Visualizamos la relación en un mapa de calor
heatmap( xtabs(~ user_age_range + distrito_station_ini, table_distrit_age), ylab="user_age_range",margins =c(10,7))#Extraemos la variable distrito inicial con movimientos que parte solo del Retiro
table_distrit_fin_age<- data %>%
dplyr::select(c(distrito_station_fin,user_age_range,distrito_station_ini)) %>%
filter(grepl("03",distrito_station_ini)==TRUE)%>%
mutate(distrito_station_fin=as.factor(distrito_station_fin),
user_age_range=as.factor(user_age_range))%>%
dplyr::select(-c(distrito_station_ini))
#Matriz de contingencia destacando frecuencia de edad (Sumatoria cols = 1)
round(prop.table(table(table_distrit_fin_age))*100, 2)## user_age_range
## distrito_station_fin 0 1 2 3 4 5 6
## 01 CENTRO 11.87 0.11 0.15 1.94 8.64 5.29 0.09
## 02 ARGANZUELA 3.66 0.09 0.07 0.44 2.32 1.76 0.01
## 03 RETIRO 14.04 0.39 0.31 2.15 8.54 8.28 0.37
## 04 SALAMANCA 6.97 0.26 0.11 0.90 5.88 3.45 0.16
## 05 CHAMARTÍN 0.95 0.05 0.01 0.15 0.87 0.58 0.00
## 06 TETUÁN 0.45 0.08 0.00 0.10 0.22 0.22 0.01
## 07 CHAMBERÍ 2.38 0.10 0.06 0.26 2.13 1.20 0.02
## 09 MONCLOA 0.79 0.00 0.03 0.08 0.59 0.45 0.00
#Realizamos también el Test de fisher
fisher.test(table(table_distrit_fin_age),simulate.p.value=TRUE)##
## Fisher's Exact Test for Count Data with simulated p-value (based
## on 2000 replicates)
##
## data: table(table_distrit_fin_age)
## p-value = 0.0004998
## alternative hypothesis: two.sided
#Visualizamos la relación
heatmap( xtabs(~ user_age_range + distrito_station_fin, table_distrit_fin_age), ylab="user_age_range",margins =c(10,7))El objetivo es detectar posibles relaciones lineales entre las variables. Si presentan un fuerte grado de correlación, buscar la forma funcional que mejor explique la variable dependiente a partir de la independiente (Análisis de regresión)
## travel_time longitude_station_ini latitude_station_ini
## Min. : 1 Min. :-3.725 Min. :40.39
## 1st Qu.: 469 1st Qu.:-3.694 1st Qu.:40.41
## Median : 684 Median :-3.688 Median :40.42
## Mean : 1250 Mean :-3.688 Mean :40.42
## 3rd Qu.: 1083 3rd Qu.:-3.678 3rd Qu.:40.42
## Max. :481814 Max. :-3.669 Max. :40.46
## longitude_station_fin latitude_station_fin distance temp_max
## Min. :-3.725 Min. :40.39 Min. : 0 Min. :30.40
## 1st Qu.:-3.694 1st Qu.:40.41 1st Qu.:1235 1st Qu.:32.70
## Median :-3.688 Median :40.42 Median :1823 Median :33.80
## Mean :-3.688 Mean :40.42 Mean :1889 Mean :34.43
## 3rd Qu.:-3.678 3rd Qu.:40.42 3rd Qu.:2397 3rd Qu.:36.70
## Max. :-3.669 Max. :40.46 Max. :5944 Max. :40.00
## temp_max_hora temp_min temp_min_hora temp_media
## Min. :15.00 Min. :16.80 Min. :5.000 Min. :23.60
## 1st Qu.:15.00 1st Qu.:19.60 1st Qu.:7.000 1st Qu.:26.10
## Median :16.00 Median :21.00 Median :8.000 Median :27.50
## Mean :15.84 Mean :21.36 Mean :7.316 Mean :27.89
## 3rd Qu.:16.00 3rd Qu.:23.00 3rd Qu.:8.000 3rd Qu.:30.10
## Max. :18.00 Max. :25.90 Max. :8.000 Max. :32.90
## viento_max viento_max_hora travel_time_mins
## Min. :11.00 Min. : 0.00 Min. : 0.017
## 1st Qu.:13.00 1st Qu.: 3.00 1st Qu.: 7.817
## Median :15.00 Median :18.00 Median : 11.400
## Mean :16.26 Mean :13.26 Mean : 20.827
## 3rd Qu.:19.00 3rd Qu.:19.00 3rd Qu.: 18.050
## Max. :25.00 Max. :23.00 Max. :8030.233
Generamos la matriz de coeficientes de correlación de Pearson y Spearman para todos los pares de variables:
# Matriz de correlacion de Pearson
corr_matrix_p <- rcorr(as.matrix(data_cont))
# Mostramos los valores de los coeficientes y no se identifica ningún valor cercano a 1 o -1
corr_matrix_p[["r"]]## travel_time longitude_station_ini
## travel_time 1.0000000000 -0.017674873
## longitude_station_ini -0.0176748728 1.000000000
## latitude_station_ini -0.0046979295 -0.019512680
## longitude_station_fin -0.0206480997 -0.097530292
## latitude_station_fin 0.0127003206 0.060140530
## distance 0.0274774723 -0.151059373
## temp_max -0.0066242061 0.006488413
## temp_max_hora 0.0093248581 -0.002565709
## temp_min -0.0064799604 0.015211328
## temp_min_hora -0.0021280376 -0.002139069
## temp_media -0.0068135284 0.011269486
## viento_max -0.0007073532 -0.002400726
## viento_max_hora -0.0019142138 0.013846186
## travel_time_mins 1.0000000000 -0.017674873
## latitude_station_ini longitude_station_fin
## travel_time -0.004697929 -0.0206480997
## longitude_station_ini -0.019512680 -0.0975302917
## latitude_station_ini 1.000000000 0.0577037050
## longitude_station_fin 0.057703705 1.0000000000
## latitude_station_fin 0.004058602 -0.0372414132
## distance 0.233569721 -0.1508242916
## temp_max 0.003365522 -0.0003444813
## temp_max_hora 0.004495453 0.0021605950
## temp_min 0.008056428 0.0043823536
## temp_min_hora 0.001616876 -0.0023939381
## temp_media 0.005966527 0.0019966595
## viento_max -0.002656509 -0.0042835918
## viento_max_hora 0.006068944 0.0094009747
## travel_time_mins -0.004697929 -0.0206480997
## latitude_station_fin distance temp_max
## travel_time 0.0127003206 0.0274774723 -0.0066242061
## longitude_station_ini 0.0601405295 -0.1510593735 0.0064884129
## latitude_station_ini 0.0040586023 0.2335697210 0.0033655215
## longitude_station_fin -0.0372414132 -0.1508242916 -0.0003444813
## latitude_station_fin 1.0000000000 0.2506925000 0.0002142818
## distance 0.2506925000 1.0000000000 -0.0065910863
## temp_max 0.0002142818 -0.0065910863 1.0000000000
## temp_max_hora 0.0051896320 -0.0100235447 0.2798889888
## temp_min 0.0063245961 0.0004258983 0.7861133440
## temp_min_hora -0.0095104846 0.0004367628 -0.1289228321
## temp_media 0.0034923128 -0.0033068938 0.9481690531
## viento_max -0.0141966883 0.0038064569 -0.2375628812
## viento_max_hora 0.0096913922 0.0176732898 -0.0548403037
## travel_time_mins 0.0127003206 0.0274774723 -0.0066242061
## temp_max_hora temp_min temp_min_hora
## travel_time 0.009324858 -0.0064799604 -0.0021280376
## longitude_station_ini -0.002565709 0.0152113281 -0.0021390692
## latitude_station_ini 0.004495453 0.0080564284 0.0016168764
## longitude_station_fin 0.002160595 0.0043823536 -0.0023939381
## latitude_station_fin 0.005189632 0.0063245961 -0.0095104846
## distance -0.010023545 0.0004258983 0.0004367628
## temp_max 0.279888989 0.7861133440 -0.1289228321
## temp_max_hora 1.000000000 0.0750279766 0.0305656189
## temp_min 0.075027977 1.0000000000 -0.1693298519
## temp_min_hora 0.030565619 -0.1693298519 1.0000000000
## temp_media 0.191772160 0.9416087152 -0.1555022294
## viento_max 0.026212462 -0.2475639319 -0.2612309556
## viento_max_hora -0.029592376 0.2267748966 0.0256348279
## travel_time_mins 0.009324858 -0.0064799604 -0.0021280376
## temp_media viento_max viento_max_hora
## travel_time -0.006813528 -0.0007073532 -0.001914214
## longitude_station_ini 0.011269486 -0.0024007262 0.013846186
## latitude_station_ini 0.005966527 -0.0026565093 0.006068944
## longitude_station_fin 0.001996659 -0.0042835918 0.009400975
## latitude_station_fin 0.003492313 -0.0141966883 0.009691392
## distance -0.003306894 0.0038064569 0.017673290
## temp_max 0.948169053 -0.2375628812 -0.054840304
## temp_max_hora 0.191772160 0.0262124623 -0.029592376
## temp_min 0.941608715 -0.2475639319 0.226774897
## temp_min_hora -0.155502229 -0.2612309556 0.025634828
## temp_media 1.000000000 -0.2556799958 0.089117332
## viento_max -0.255679996 1.0000000000 0.040719742
## viento_max_hora 0.089117332 0.0407197419 1.000000000
## travel_time_mins -0.006813528 -0.0007073532 -0.001914214
## travel_time_mins
## travel_time 1.0000000000
## longitude_station_ini -0.0176748728
## latitude_station_ini -0.0046979295
## longitude_station_fin -0.0206480997
## latitude_station_fin 0.0127003206
## distance 0.0274774723
## temp_max -0.0066242061
## temp_max_hora 0.0093248581
## temp_min -0.0064799604
## temp_min_hora -0.0021280376
## temp_media -0.0068135284
## viento_max -0.0007073532
## viento_max_hora -0.0019142138
## travel_time_mins 1.0000000000
# Heatmap Pearson
corrplot(corr_matrix_p$r, type = "upper", order = "alphabet", tl.col = "black", win.asp = 0.8, tl.cex = 0.8 , title="Heatmap - Pearson", mar=c(0,0,1,0))# Matriz de correlación de Spearman
corr_matrix_s <- rcorr(as.matrix(data_cont), type="spearman")
# Mostramos los valores de los coeficientes y no se identifica ningún valor cercano a 1 o -1
corr_matrix_s[["r"]]## travel_time longitude_station_ini
## travel_time 1.000000000 -0.094961698
## longitude_station_ini -0.094961698 1.000000000
## latitude_station_ini 0.035932870 0.024231486
## longitude_station_fin -0.102620135 -0.073095437
## latitude_station_fin 0.106031785 0.040563783
## distance 0.469928266 -0.117000020
## temp_max -0.032810472 0.003177315
## temp_max_hora 0.006815319 -0.001746893
## temp_min -0.040524887 0.011266963
## temp_min_hora -0.006415498 0.001406170
## temp_media -0.039473714 0.008964384
## viento_max 0.011486869 -0.001551442
## viento_max_hora -0.025575981 0.009992054
## travel_time_mins 1.000000000 -0.094961698
## latitude_station_ini longitude_station_fin
## travel_time 0.035932870 -0.1026201352
## longitude_station_ini 0.024231486 -0.0730954373
## latitude_station_ini 1.000000000 0.0381286270
## longitude_station_fin 0.038128627 1.0000000000
## latitude_station_fin 0.018471512 0.0039950635
## distance 0.095878553 -0.1183316867
## temp_max 0.006916440 -0.0004112203
## temp_max_hora 0.004927458 0.0026335360
## temp_min 0.010850878 0.0020978586
## temp_min_hora 0.005432258 -0.0001942805
## temp_media 0.010821959 0.0020264022
## viento_max -0.006910896 -0.0020979001
## viento_max_hora 0.004396632 0.0078291242
## travel_time_mins 0.035932870 -0.1026201352
## latitude_station_fin distance temp_max
## travel_time 0.106031785 0.4699282661 -0.0328104722
## longitude_station_ini 0.040563783 -0.1170000197 0.0031773152
## latitude_station_ini 0.018471512 0.0958785529 0.0069164404
## longitude_station_fin 0.003995064 -0.1183316867 -0.0004112203
## latitude_station_fin 1.000000000 0.1145705635 0.0079022528
## distance 0.114570564 1.0000000000 -0.0032922756
## temp_max 0.007902253 -0.0032922756 1.0000000000
## temp_max_hora 0.007735662 -0.0122418514 0.1528658527
## temp_min 0.009723055 -0.0010792157 0.7655483162
## temp_min_hora -0.003054865 0.0008732965 -0.0461419882
## temp_media 0.009662459 -0.0038717199 0.9303800647
## viento_max -0.009691485 0.0040489509 -0.2393712274
## viento_max_hora 0.011372480 0.0137994236 -0.0805360277
## travel_time_mins 0.106031785 0.4699282661 -0.0328104722
## temp_max_hora temp_min temp_min_hora
## travel_time 0.006815319 -0.040524887 -0.0064154981
## longitude_station_ini -0.001746893 0.011266963 0.0014061699
## latitude_station_ini 0.004927458 0.010850878 0.0054322583
## longitude_station_fin 0.002633536 0.002097859 -0.0001942805
## latitude_station_fin 0.007735662 0.009723055 -0.0030548652
## distance -0.012241851 -0.001079216 0.0008732965
## temp_max 0.152865853 0.765548316 -0.0461419882
## temp_max_hora 1.000000000 0.106337542 0.1566595557
## temp_min 0.106337542 1.000000000 0.0112661711
## temp_min_hora 0.156659556 0.011266171 1.0000000000
## temp_media 0.139980278 0.931195119 -0.0059619988
## viento_max 0.025417408 -0.218518385 -0.2680092973
## viento_max_hora 0.111698285 0.088291250 0.0274262549
## travel_time_mins 0.006815319 -0.040524887 -0.0064154981
## temp_media viento_max viento_max_hora
## travel_time -0.039473714 0.011486869 -0.025575981
## longitude_station_ini 0.008964384 -0.001551442 0.009992054
## latitude_station_ini 0.010821959 -0.006910896 0.004396632
## longitude_station_fin 0.002026402 -0.002097900 0.007829124
## latitude_station_fin 0.009662459 -0.009691485 0.011372480
## distance -0.003871720 0.004048951 0.013799424
## temp_max 0.930380065 -0.239371227 -0.080536028
## temp_max_hora 0.139980278 0.025417408 0.111698285
## temp_min 0.931195119 -0.218518385 0.088291250
## temp_min_hora -0.005961999 -0.268009297 0.027426255
## temp_media 1.000000000 -0.226439217 -0.012671728
## viento_max -0.226439217 1.000000000 0.176807966
## viento_max_hora -0.012671728 0.176807966 1.000000000
## travel_time_mins -0.039473714 0.011486869 -0.025575981
## travel_time_mins
## travel_time 1.000000000
## longitude_station_ini -0.094961698
## latitude_station_ini 0.035932870
## longitude_station_fin -0.102620135
## latitude_station_fin 0.106031785
## distance 0.469928266
## temp_max -0.032810472
## temp_max_hora 0.006815319
## temp_min -0.040524887
## temp_min_hora -0.006415498
## temp_media -0.039473714
## viento_max 0.011486869
## viento_max_hora -0.025575981
## travel_time_mins 1.000000000
# Heatmap Spearman
corrplot(corr_matrix_s$r, type = "upper", order = "alphabet", tl.col = "black", win.asp = 0.8, tl.cex = 0.8, title="Heatmap - Spearman", mar=c(0,0,1,0))ggplot( data_cont , aes(x = distance, y = travel_time_mins)) +
geom_point() +
theme(plot.subtitle = element_text(vjust = 1),
plot.caption = element_text(vjust = 1),
panel.grid.major = element_line(colour = "orange"),
axis.title = element_text(face = "bold"),
axis.text = element_text(face = "bold"),
plot.title = element_text(face = "bold", colour = "chocolate")) +
labs(title = "Estudio Correlación Distancia - Tiempo Viaje", subtitle = "")Eliminamos Outliers que como vimos antes son duraciones superiores a 1,30 horas segundos y descartando a los trabajadores
y <- data %>%
select(distance, travel_time_mins,user_type) %>%
filter(travel_time_mins < 87,
distance > 0,
user_type=='1'|user_type=='2') %>%
dplyr::select(-c(user_type))
ggplot( y , aes(x = distance, y = travel_time_mins)) +
geom_point() +
geom_smooth(method=lm) +
theme(plot.subtitle = element_text(vjust = 1),
plot.caption = element_text(vjust = 1),
panel.grid.major = element_line(colour = "orange"),
axis.title = element_text(face = "bold"),
axis.text = element_text(face = "bold"),
plot.title = element_text(face = "bold", colour = "chocolate")) +
labs(title = "Estudio Correlación Distancia - Tiempo Viaje", subtitle = "")#Extraemos la variable distrito inicial eliminado Retiro
table_travel_user <- data %>%
dplyr::select(c(travel_time_mins, user_type)) %>%
mutate(user_type = as.numeric(user_type))
#calculamos coeficiente de pearson
cor.test(table_travel_user$travel_time_mins, table_travel_user$user_type, method = c("pearson"))##
## Pearson's product-moment correlation
##
## data: table_travel_user$travel_time_mins and table_travel_user$user_type
## t = 44.664, df = 60437, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1710243 0.1864597
## sample estimates:
## cor
## 0.178753
#filtramos para los valores 1 y 2, eliminando trayectos de la empresa
z <- table_travel_user %>%
filter(user_type=='1'|user_type=='2')
cor.test(z$travel_time_mins, z$user_type, method = c("pearson"))##
## Pearson's product-moment correlation
##
## data: z$travel_time_mins and z$user_type
## t = 86.378, df = 54738, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3389511 0.3536958
## sample estimates:
## cor
## 0.3463448
No hay correlación entre el tipo de usuario y el tiempo. Se puede ver que eliminado los trayectos de los trabajadores, subre ligeramente el coeficiente de pearson, ya que los trayecto son más duraderos.
#Extraemos la variable distrito inicial eliminado Retiro
table_travel_age <- data %>%
dplyr::select(c(travel_time_mins, user_age_range)) %>%
mutate(user_age_range = as.numeric(user_age_range))
#calculamos coeficiente de pearson
cor.test(table_travel_age$user_age_range, table_travel_age$travel_time_mins, method = c("pearson"))##
## Pearson's product-moment correlation
##
## data: table_travel_age$user_age_range and table_travel_age$travel_time_mins
## t = -5.1455, df = 60437, p-value = 2.676e-07
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.02889337 -0.01295548
## sample estimates:
## cor
## -0.02092575
Al filtrar los usuarios de edad desconocidos el p_value pasa de negativo a ser >0,5, aunque el test no refleja correlación directa, si que parece que a mayor edad menor tiempo
#filtramos los desconocidos
table_trave_age_filter <- table_travel_age%>%
filter(user_age_range!=0)
cor.test(table_trave_age_filter$user_age_range, table_trave_age_filter$travel_time_mins, method = c("pearson"))##
## Pearson's product-moment correlation
##
## data: table_trave_age_filter$user_age_range and table_trave_age_filter$travel_time_mins
## t = 0.51208, df = 34963, p-value = 0.6086
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.007743348 0.013220007
## sample estimates:
## cor
## 0.00273863
Estudio por hora de la posible correlación entre la hora de temperatura máxima e inicio del trayecto.
#Generamos la tabla con dos variables
table_hour_max <- data %>%
dplyr::select(c(unplug_hour, temp_max_hora)) %>%
mutate (unplug_hour= as.numeric(unplug_hour),
temp_max_hora = as.numeric(temp_max_hora))
#Realizamos también el Test de Pearson
cor.test(table_hour_max$unplug_hour, table_hour_max$temp_max_hora, method = c("pearson"))##
## Pearson's product-moment correlation
##
## data: table_hour_max$unplug_hour and table_hour_max$temp_max_hora
## t = -4.5098, df = 60437, p-value = 6.501e-06
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.02630995 -0.01037044
## sample estimates:
## cor
## -0.01834136
Se utiliza PAM (Partición Alrededor de Medoids) que es un algoritmo de clasificación del tipo k-medoids. Escoge datapoints como centros y trabaja con una métrica arbitraria de distancias entre datapoints.
vars_cluster <- data %>%
filter(between (unplug_date,'01/08/2018','28/08/2018')) %>%
filter( user_type == '1') %>%
select(idunplug_station_ini,cp_ini,idplug_station_fin,cp_fin,travel_time,user_type,
user_age_range,unplug_hour,distance,day, bank_holiday_flg) %>%
mutate(
day = case_when(
day =='L' ~ 0,
day =='M' ~ 1,
day =='X' ~ 2,
day =='J' ~ 3,
day =='V' ~ 4,
day =='S' ~ 5,
day =='D' ~ 6),
user_age_range=as.factor(user_age_range),
unplug_hour=as.factor(unplug_hour),
user_type=as.factor(user_type),
idunplug_station_ini=as.factor(idunplug_station_ini),
idplug_station_fin=as.factor(idplug_station_fin),
cp_ini=as.factor(cp_ini),
cp_fin=as.factor(cp_fin),
bank_holiday_flg=as.factor(bank_holiday_flg))
#filtramos para quedarnos con los tipos 1 y 2
#Selecionamos variables
data_cluster <- vars_cluster %>%
select(idunplug_station_ini,cp_ini,idplug_station_fin,cp_fin,travel_time,user_type,
user_age_range,unplug_hour,distance,day,bank_holiday_flg)
rm (vars_cluster)En esta ejecucion del algoritmo se utilizan los datos con outliers (sin filtrado previo). Hacemos un subconjunto de 2500 muestras y calculamos el número de clusters recomendados para el Kmeans con el método del codo.
#Subsetting para clasterizar.
data_cluster_sample <- sample_n(data_cluster, 2500)
#Num cluster - Elbow method for KMEANS
fviz_nbclust(x = data_cluster_sample, FUNcluster = kmeans, method = "wss", k.max = 20) +
geom_vline(xintercept = 5, linetype = 2)Calculamos la distancia de Gower más representativa para el algoritmo PAM y utilizamos el método Silhouette, donde se escoge el primer pico cómo número posible de clusters.
data_cluster <- data_cluster %>%
select(idunplug_station_ini,cp_ini,idplug_station_fin,cp_fin,travel_time,user_type,
user_age_range,unplug_hour,distance,day)
rm (vars_cluster)## Warning in rm(vars_cluster): object 'vars_cluster' not found
# PAM - Gower Distance
gower_dist <- daisy(data_cluster_sample, metric = "gower")
#2 clusters has the highest silhouette width. Let’s pick k = 5
#search for a number of clusters
sil_width <- c(NA)
for(i in 2:8){
pam_fit <- pam(gower_dist, diss = TRUE, k = i)
sil_width[i] <- pam_fit$silinfo$avg.width
}
plot(1:8, sil_width,
xlab = "Number of clusters",
ylab = "Silhouette Width")
lines(1:8, sil_width)Ejecutamos el algoritmo PAM para 4 clusteres
#executing the algorithm with 4 clusters
k <- 4
pam_fit <- pam(gower_dist, diss = TRUE, k)
pam_results <- data_cluster_sample %>%
mutate(cluster = pam_fit$clustering) %>%
group_by(cluster) %>%
do(the_summary = summary(.))
#Visualization in a lower dimensional space
tsne_obj <- Rtsne(gower_dist, is_distance = TRUE)
tsne_data <- tsne_obj$Y %>%
data.frame() %>%
setNames(c("X", "Y")) %>%
mutate(cluster = factor(pam_fit$clustering))
ggplot(aes(x = X, y = Y), data = tsne_data) +
geom_point(aes(color = cluster))Se observa que la clasterización no es muy buena con 4 clusters. Probamos con 5 que es lo que nos recomiendan ambos métodos de distancia.
Se obtiene una mejor representación, aunque debido a la muestra de 2500 muestras a veces la clasificación no es muy adecuada..
#The results looks like 3-4cluster, so we enter again.
k <- 5
pam_fit <- pam(gower_dist, diss = TRUE, k)
pam_results <- data_cluster_sample %>%
mutate(cluster = pam_fit$clustering) %>%
group_by(cluster) %>%
do(the_summary = summary(.))Visualización del modelo con 5 clusters.
#Visualization in a lower dimensional space
tsne_obj <- Rtsne(gower_dist, is_distance = TRUE)
tsne_data <- tsne_obj$Y %>%
data.frame() %>%
setNames(c("X", "Y")) %>%
mutate(cluster = factor(pam_fit$clustering))
ggplot(aes(x = X, y = Y), data = tsne_data) +
geom_point(aes(color = cluster))Resumen de resultados los 5 clusters.
## [[1]]
## idunplug_station_ini cp_ini idplug_station_fin cp_fin
## 80 : 30 28007 :102 64 : 96 28001 :108
## 74 : 27 28009 : 62 91 : 17 28014 : 55
## 86 : 19 28014 : 60 73 : 16 28009 : 54
## 84 : 17 28012 : 55 100 : 15 28007 : 53
## 73 : 16 28004 : 35 80 : 14 28012 : 49
## 102 : 11 28006 : 24 74 : 13 28006 : 31
## (Other):317 (Other): 99 (Other):266 (Other): 87
## travel_time user_type user_age_range unplug_hour distance
## Min. : 157 1:437 0:281 18 : 46 Min. : 0
## 1st Qu.: 477 1: 0 19 : 35 1st Qu.:1231
## Median : 664 2: 3 20 : 34 Median :1894
## Mean : 1038 3: 37 11 : 26 Mean :1883
## 3rd Qu.: 1026 4: 71 17 : 26 3rd Qu.:2391
## Max. :20611 5: 43 22 : 24 Max. :5661
## 6: 2 (Other):246
## day bank_holiday_flg cluster
## Min. :0.00 0: 77 Min. :1
## 1st Qu.:5.00 1:360 1st Qu.:1
## Median :5.00 Median :1
## Mean :4.76 Mean :1
## 3rd Qu.:6.00 3rd Qu.:1
## Max. :6.00 Max. :1
##
##
## [[2]]
## idunplug_station_ini cp_ini idplug_station_fin cp_fin
## 90 :118 28009 :346 99 : 51 28004 :138
## 65 : 46 28014 : 29 73 : 14 28012 : 55
## 75 : 40 28004 : 28 86 : 14 28014 : 46
## 66 : 35 28006 : 19 9 : 12 28006 : 42
## 102 : 29 28012 : 13 170 : 11 28009 : 40
## 101 : 27 28045 : 7 31 : 11 28007 : 31
## (Other):177 (Other): 30 (Other):359 (Other):120
## travel_time user_type user_age_range unplug_hour distance
## Min. : 116.0 1:472 0:119 15 : 81 Min. : 0
## 1st Qu.: 387.2 1: 1 9 : 45 1st Qu.:1081
## Median : 557.5 2: 1 8 : 34 Median :1694
## Mean : 761.8 3: 27 14 : 31 Mean :1748
## 3rd Qu.: 835.8 4:242 20 : 30 3rd Qu.:2281
## Max. :7568.0 5: 80 18 : 28 Max. :5141
## 6: 2 (Other):223
## day bank_holiday_flg cluster
## Min. :0.000 0:427 Min. :2
## 1st Qu.:1.000 1: 45 1st Qu.:2
## Median :2.000 Median :2
## Mean :2.176 Mean :2
## 3rd Qu.:3.000 3rd Qu.:2
## Max. :6.000 Max. :2
##
##
## [[3]]
## idunplug_station_ini cp_ini idplug_station_fin cp_fin
## 79 : 76 28007 :236 90 :111 28009 :309
## 74 : 49 28012 : 67 65 : 39 28012 : 52
## 80 : 34 28014 : 48 66 : 38 28007 : 48
## 78 : 32 28004 : 45 75 : 36 28006 : 41
## 73 : 19 28006 : 45 102 : 25 28014 : 37
## 86 : 16 28045 : 21 107 : 17 28001 : 31
## (Other):380 (Other):144 (Other):340 (Other): 88
## travel_time user_type user_age_range unplug_hour distance
## Min. : 1.0 1:606 0:128 17 : 63 Min. : 0
## 1st Qu.: 494.5 1: 0 8 : 60 1st Qu.:1502
## Median : 667.5 2: 3 7 : 48 Median :2013
## Mean : 868.6 3: 41 19 : 40 Mean :2090
## 3rd Qu.: 893.5 4:339 20 : 38 3rd Qu.:2533
## Max. :18543.0 5: 93 21 : 37 Max. :5779
## 6: 2 (Other):320
## day bank_holiday_flg cluster
## Min. :0.000 0:564 Min. :3
## 1st Qu.:1.000 1: 42 1st Qu.:3
## Median :2.000 Median :3
## Mean :2.076 Mean :3
## 3rd Qu.:3.000 3rd Qu.:3
## Max. :6.000 Max. :3
##
##
## [[4]]
## idunplug_station_ini cp_ini idplug_station_fin cp_fin
## 64 :147 28001 :172 84 : 80 28014 :236
## 73 : 19 28012 : 60 91 : 43 28009 : 59
## 99 : 19 28009 : 59 85 : 28 28007 : 47
## 65 : 15 28006 : 42 73 : 27 28006 : 39
## 74 : 15 28014 : 39 69 : 21 28012 : 31
## 86 : 13 28004 : 36 72 : 18 28016 : 28
## (Other):326 (Other):146 (Other):337 (Other):114
## travel_time user_type user_age_range unplug_hour distance
## Min. : 113.0 1:554 0:369 13 : 51 Min. : 0
## 1st Qu.: 428.0 1: 0 19 : 45 1st Qu.:1172
## Median : 610.0 2: 3 20 : 39 Median :1693
## Mean : 891.0 3: 26 15 : 35 Mean :1770
## 3rd Qu.: 907.8 4: 92 7 : 35 3rd Qu.:2217
## Max. :17249.0 5: 60 18 : 32 Max. :5618
## 6: 4 (Other):317
## day bank_holiday_flg cluster
## Min. :0.000 0:497 Min. :4
## 1st Qu.:2.000 1: 57 1st Qu.:4
## Median :3.000 Median :4
## Mean :2.684 Mean :4
## 3rd Qu.:4.000 3rd Qu.:4
## Max. :6.000 Max. :4
##
##
## [[5]]
## idunplug_station_ini cp_ini idplug_station_fin cp_fin
## 91 : 57 28014 :162 74 : 86 28007 :212
## 84 : 32 28012 : 42 79 : 32 28009 : 36
## 69 : 26 28009 : 41 78 : 28 28012 : 33
## 99 : 18 28004 : 38 80 : 24 28006 : 31
## 85 : 16 28007 : 32 100 : 11 28014 : 18
## 74 : 12 28006 : 30 73 : 11 28001 : 16
## (Other):270 (Other): 86 (Other):239 (Other): 85
## travel_time user_type user_age_range unplug_hour distance
## Min. : 115.0 1:431 0:105 14 : 75 Min. : 0
## 1st Qu.: 466.5 1: 3 15 : 31 1st Qu.:1299
## Median : 654.0 2: 4 19 : 29 Median :1862
## Mean : 872.2 3: 18 20 : 29 Mean :1976
## 3rd Qu.: 948.0 4: 64 8 : 27 3rd Qu.:2420
## Max. :10479.0 5:233 16 : 26 Max. :5552
## 6: 4 (Other):214
## day bank_holiday_flg cluster
## Min. :0.000 0:386 Min. :5
## 1st Qu.:0.000 1: 45 1st Qu.:5
## Median :1.000 Median :5
## Mean :1.838 Mean :5
## 3rd Qu.:3.000 3rd Qu.:5
## Max. :6.000 Max. :5
##
La parte positiva del modelo de PAM es que obtenemos unos resultados que aunque habría que ejecutar con una muestra más amplia o mejor definida, nos da unos patrones que nos permitirían explicar de manera clara a negocio los diferentes clusters.
Se utiliza CLARA (Clustering Large Applications) es un algoritmo de clasificacion basado en K-Means que utiliza tecnicaas de muestreo para caracterizar grandes datasets. Toma un numero pre-determinado de clusters.
# Seleccion de variables para el modelo
vars_modelo <- c('idunplug_station_ini','cp_ini','idplug_station_fin','cp_fin','travel_time','user_type', 'user_age_range','unplug_hour','day','distance')En esta ejecucion del algoritmo se utilizan los datos con outliers (sin filtrado previo) y 6 clusters.
# CON OUTLIERS
# dataset
data_modelo <- data %>%
select(vars_modelo)
# Convertir data types
data_modelo[sapply(data_modelo, is.character)] <- lapply(data_modelo[sapply(data_modelo, is.character)], as.factor)
# Compute CLARA
# clara_result <- clara(data_modelo, 10, samples = 2500, pamLike = TRUE)
load(file = 'clara_result.RData')
# Add the point classifications to the original data, use this:
data_modelo_clusters <- cbind(data_modelo, cluster = clara_result$cluster)
head(data_modelo_clusters, n = 5)## idunplug_station_ini cp_ini idplug_station_fin cp_fin travel_time
## 1 99 28004 103 28009 291
## 2 86 28012 77 28009 549
## 3 100 28006 52 28012 632
## 4 44 28012 102 28009 847
## 5 165 28046 91 28014 1289
## user_type user_age_range unplug_hour day distance cluster
## 1 1 5 1 J 1334.102 1
## 2 1 4 1 J 1970.717 2
## 3 1 0 1 J 2293.919 3
## 4 1 0 1 J 2319.953 3
## 5 1 0 1 J 1902.985 3
En el grafico se observa que el cluster-1 agrupa todos los aoutliers del dataset.
En esta ejecucion del algortimo se toma el datatset original filtrado por quartiles extremos con la funcion boxplot.stats()$out. Se toman 5 clusters.
# SIN OUTLIERS
# datatset
data_modelo_nooutliers <- data_no_outliers %>%
select(vars_modelo)
# Convertir data types
data_modelo_nooutliers[sapply(data_modelo_nooutliers, is.character)] <- lapply(data_modelo_nooutliers[sapply(data_modelo_nooutliers, is.character)], as.factor)
# Compute CLARA
# clara_result_nooutliers <- clara(data_modelo_nooutliers, 10, samples = 2500, pamLike = TRUE)
load(file = 'clara_result_nooutliers.RData')
# Add the point classifications to the original data, use this:
data_modelo_nooutliers_clusters <- cbind(data_modelo_nooutliers, cluster = clara_result_nooutliers$cluster)
head(data_modelo_nooutliers_clusters, n = 4)## idunplug_station_ini cp_ini idplug_station_fin cp_fin travel_time
## 1 99 28004 103 28009 291
## 2 86 28012 77 28009 549
## 3 100 28006 52 28012 632
## 4 44 28012 102 28009 847
## user_type user_age_range unplug_hour day distance cluster
## 1 1 5 1 J 1334.102 1
## 2 1 4 1 J 1970.717 2
## 3 1 0 1 J 2293.919 2
## 4 1 0 1 J 2319.953 3
Se observa que los grupos est’an agrupados entre si, pero no hay mucha distancia entre clusters.
-Los picos de uso de Bicimad se observan muy marcados con el inicio y fin de la jornada laboral, podría ser interesantes realizar acciones comerciales dirigidas hacia esas horas
-Excluyendo los de edad desconocida, el 86% de usuarios son adultos (27-65 años)
-Uso constante casi rutinario a nivel de semana, cayendo en los fines de semana, más pronunciado esto en los usuarios de edad adulta (27-65 años)
-En lo referente a los clusters, aunque no se observan grandes diferencias entre los centroides, si aparecen bien definidos y explicados por las variables seleccionadas
Comparando con los datos de Enero: - El uso del servicio baja en las horas de mayor calor 15-18 horas, se intensifica entre las 7-9 de la mañana. - En verano el uso entre las 19-22 horas se amplía - Los trayectos son más cortos, destacan muchos movimiento en el propio retiro o con estaciones centrales.
PROXIMOS PASOS
-Utilizar herramientas desktop para no tener problemáticas de memoria
-Cruzar con datasets de eventos ligado a zonas
-Volver a ejecutar el proyecto con datos de otros periodos, pues están algo sesgados por la temporalidad, para poder ver como se adaptan.